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We present a new linearization of T-Matrix and Mie computations for light scattering 
by non-spherical and spherical particles, respectively. In addition to the usual extinc- 
tion and scattering cross-sections and the scattering matrix outputs, the linearized 
models will generate analytical derivatives of these optical properties with respect to 
the real and imaginary parts of the particle refractive index, and (for non-spherical 
scatterers) with respect to the “shape” parameter (the spheroid aspect ratio, cylinder 
diameter/height ratio, Chebyshev particle deformation factor). These derivatives are 
based on the essential linearity of Maxwell’s theory. Analytical derivatives are also 
available for polydisperse particle size distribution parameters such as the mode radius. 
The T-matrix formulation is based on the NASA Goddard Institute for Space Studies 
FORTRAN 77 code developed in the 1990s. The linearized scattering codes presented 
here are in FORTRAN 90 and will be made publicly available. 

© 201 1 Elsevier Ltd. Ah rights reserved. 


1. Introduction 

The generation of accurate light scattering properties 
for spherical and non-spherical particles is extremely 
important for many applications in a wide variety of 
physical science disciplines. Of particular importance are 
methods based on direct numerical solutions of Maxwell’s 
equations of electrodynamics. The first accurate light 
scattering calculations for spherical particles date back 
to the pioneering work of Mie and Lorenz (see [1-3] 
for reviews and perspective). There are many Mie codes 
available in the public domain; in this work, our basis is a 
model generated in the 1980s by a Dutch group [4]. 

For non-spherical particles, there are several methods 
for computing optical properties; for a review, see [5]. Of 
these methods, the T-matrix approach first conceived by 
Waterman [6] has been developed extensively in the 
last two decades for a huge variety of applications; the 
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data-base review [7] is useful in this regard. In this work, 
our starting point is the popular and widely available 
T-matrix code disseminated by the NASA Goddard Institute 
for Space Studies (GISS) group [8,9]. Other codes are 
reviewed in [10]. The NASA-GISS code is applicable to 
randomly oriented spheroids, circular cylinders and Cheby- 
shev particles. The reader is referred to two papers for 
details: Ref. [8] presents a review of the theory, while Ref. 
[9] presents a description of the FORTRAN 77 code for 
numerical computations. 

With changing climate dynamics, it has become 
important to obtain accurate quantitative information 
on aerosol optical properties on a global scale [11] from 
both dedicated ground-based and space-borne instru- 
ments [12]. To date, retrievals of aerosol optical thickness 
are commonplace for many remote sensors. But in the 
absence of polarimetric measurements, it is difficult 
to obtain additional information (such as aerosol single 
scattering albedo) that is important for estimates of 
aerosol climate forcing. With the recent deployment of 
polarimetric sensors such as the Research Scanning 
Polarimeter (RSP) [13], the potential for extending and 
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improving the retrieval of aerosol parameters to include 
absorption properties was clearly demonstrated in a number 
of studies (see for example [14,15]). The recent tragedy of 
the aborted GLORY Mission [16] has deprived the com- 
munity of a valuable tool for space-borne aerosol detec- 
tion. However the RSP instrument will continue to be 
deployed from air-borne platforms [17]. 

A 2005 study on GOME-2 measurements [18] demon- 
strated the feasibility of deriving microphysical aerosol 
parameters (refractive index, size distribution parameters) 
in addition to the more usual macrophysical optical 
properties such as aerosol extinction and scattering pro- 
files. The retrieval was based on a forward model com- 
prising a linearized vector radiative transfer (RT) code 
acting alongside a linearized Mie model. The latter gen- 
erates analytic partial derivatives of optical properties 
with respect to microphysical aerosol parameters. This 
type of combination tool is particularly useful for inverse 
and sensitivity algorithms requiring analytic Jacobians (of 
atmospheric parameters) in addition to the usual radia- 
tion field simulations. 

Another study for the OCO instrument used a similar 
approach [19], this time in connection with the retrieval 
of XCO 2 columns from the weak and strong CO 2 bands 
(1.60 and 2.04 pm, OCO also samples the O 2 A band); 
aerosol characterization is an essential part of this retrie- 
val, given the requirement to obtain CO 2 estimates at 
l-3ppmv accuracy [20]. Rather than overburden the 
forward model by specifying macrophysical aerosol opti- 
cal properties in every layer, this study used a parameter- 
ized tropospheric aerosol formulation with simple expo- 
nential, linear or Gaussian loading profiles, and a handful 
of microphysical Mie-based aerosol properties. The latter 
are then retrieved along with the total loading and 
another parameter (such as the exponential relaxation 
constant) characterizing the loading profile. This method 
allows for a better characterization of aerosol uncertainty 
as a source of forward model error in the CO 2 retrieval. 

However, at wavelengths in and around the O 2 A 
absorption band, the vertical profile of aerosol scattering 
is important for the accurate simulation of radiance and 
polarization at top-of-atmosphere [21]. By analogy with 
UV aerosol retrieval algorithms in which the profile of 
Rayleigh scattering is used to calibrate the profile of 
(high-altitude) absorbing aerosols [22], Zeng et al. [21] 
showed that, for polarization measurements at the O 2 A 
band, the profile of O 2 absorption may be used to calibrate 
the profile of scattering aerosols [21]. 

The increasing need for knowledge of 3D aerosol 
optical properties for both climate studies and satellite 
remote sensing applications requires the development 
of accurate measurements of all Stokes parameters for 
characterizing aerosol scattering, as well as the develop- 
ment of modeling tools that can rapidly and accurately 
simulate the sensitivity of the four Stokes parameters of 
the scattered light to changes in aerosol microphysical 
parameters. In line with this goal, the present authors 
have constructed a general tool for aerosol property 
retrieval based on the linearized VLIDORT polarization 
RT model [23] and the linearized Mie code outlined in 
this paper. 


It is well known that non-spherical dust particles are 
omnipresent in the atmosphere, and these have different 
phase functions compared to those for spherical particles [8]; 
such differences can lead to significant errors in ground- 
based or satellite-based retrieval of aerosol optical thick- 
ness and other aerosol parameters, as demonstrated by 
Refs. [24,25], and references therein. Hence, the sensitiv- 
ity of Stokes parameters to changes of (non-spherical) 
particle characteristics is important, and this sensitivity 
can be provided by the combination of VLIDORT and the 
linearized T-matrix code for the remote sensing of aerosol 
properties. 

The Mie code was linearized independently in [18,26] 
as well as by one of the present authors [R. Spurr, 2004, 
unpublished note]. Here we present a new linearization of 
the T-matrix formulation. For individual particles we 
show that the T-matrix theory is analytically differenti- 
able with respect to the three microphysical variables — 
the real and imaginary parts and mi of the particle 
refractive index mc=mr-\-imi, and the particle deforma- 
tion characteristic or shape parameter s (for spheroids, 
this is the ratio of the semi-axes; for cylinders, the 
diameter/height ratio; for Chebyshev particles, the defor- 
mation parameter). 

In Section 2 we present an overview of the T-matrix 
formulation, including a definition of the linearization 
process for the T-matrix itself. In Section 3 we discuss 
in detail analytic differentiation of the vector spherical 
functions and integrals over the particle surface areas 
with respect to mr, mi and s. Section 4 deals with poly- 
disperse linearizations with respect to parameters char- 
acterizing equivalent-sphere particle size distribution. In 
Section 5, we present some results for extinction and 
scattering cross-sections and scattering matrices and their 
linearizations. Section 6 gives a brief digest of the new 
Fortran 90 computer code for this linearization. 

2. Basic definitions and the linearization principle 

2.1. Optical properties and linearizations 

We consider the scattering of light by spherical parti- 
cles (Mie) or non-spherical particles with an axis of 
rotational symmetry. Particles are assumed to be ran- 
domly oriented and to scatter independently. The scatter- 
ing is characterized by the extinction cross-section per 
particle Cext, the scattering cross-section Csca per particle, 
and the 4x4 normalized scattering matrix F(0) for 
scattering angle 0 [27]. These quantities are ensemble- 
averaged over all orientations. The absorption cross-sec- 
tion is Cabs =Cext- Csca. and the single scattering albedo is 
^ — Cscal Cext‘ 

In the conventional phenomenological description of 
far-held scattering by a volume element dv, the scattering 
and incident Stokes 4-vectors hca and Imc are related 
through 

ry ^sca^O dt'F(0)/jf|c* (1) 

where R is the distance to a far-held observation point, 
and no the particle number density. As noted in recent 
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work by Mishchenko (see for example [28]), expressions 
such as (1) are properly valid when certain well-defined 
conditions are observed, and for this reason we work only 
with optical properties Cext, Csca and F(0). 

For the particles considered in this paper, F(0) has the 
form 


F(0) = 


/ai(0) 

bi(0) 

0 

0 


bi(0) 

02 ( 0 ) 

0 

0 


0 0 \ 

0 0 

03(0) b2(0) 

-f)2(0) 04(0 ) ) 


( 2 ) 


where there are only six independent quantities (four for 
Mie scattering). It is convenient (and more efficient) for 
most applications to use expansions of these F-matrix 
entries in terms of generalized spherical functions P\nn(^) 

LM LM 

Oi (0) = Y) ^oo(cos 0); 04(0) =Y «4Poo(cos 0); 

1 = 0 1 = 0 

(3a) 

LM 

02(0) ± 03 ( 0 ) = + « 3 )Pk ±2(cos0); (3b) 

1 = 0 


LM LM 

bi(0)= 5 ^MPo2(cos 0); b2(0)= X^i?2Po2(cos0). (3c) 

/ = 0 1 = 0 

The (1, 1) entry is the phase function, represented as 
an expansion in terms of ordinary Legendre polynomials; 
it is normalized to unity. We note also the asymmetry 
parameter: g = 1 /3a]. For more details, see for example [27]. 

The basic set of optical properties for a single 
particle is then 

W = {Cexf*Csca*<^i »^2i- ("^) 

For polydisperse ensembles, we must average over the 
particle size distribution (PSD). If n{r,v)dr is the number of 
particles in the range [r, r+dr], n and T 2 are the minimum 
and maximum such radii, and N(v) is the particle number 
density, then the polydisperse cross-sections < Cext > » < Csca > 
and the expansion coefficient sets <y,> (where {yd is one of 
{a^,a^,a^,aJj,jS^,jS 2 }) are given by 

<^ext(r)nir,v)dr; Csca = j^^ f Csca(r)n(r,v)dr; 

(5a) 

yi=d~ f 7 i(r)Csca(r)n(r.v)dr. (5b) 

Csca Jn 

Flere, vector v is shorthand for the set of parameters 
characterizing the PSD; for example v={rg, Sg] for a lognor- 
mal distribution with mode radius Tg and standard deviation 
Sg. Integrations are done numerically, usually with Gauss- 
Legendre quadrature. 

Bimodal distributions are common in aerosol retrie- 
vals; in this case we have separate sets and of 
monodisperse optical properties, plus associated PSDs 
n^^\r) and Total polydisperse cross-sections and 

expansion coefficients are given by 

Cext, sea = fC^ext,sca + (^ ~/)^S,sca ’ 


/cWy<n(l-/)Cg>yp) 

/cS+(i-/)cg> 


(6b) 


Flere, / = is the fractional number density 

corresponding to PSD Often, the two distributions 

are of the same form (e.g. both lognormal), and sometimes 
PSD properties will be shared, e.g. a common lognormal 
standard deviation but different mode radii [18]. 

A linearized T-matrix or Mie scattering model will not 
only produce the above set of properties in Eqs. (4), (5a) 
and (5b), but also their analytic partial derivatives (i) with 
respect to the individual-particle microphysical proper- 
ties nir, rrii and £, and (ii) with respect to any member Vt of 
the set of parameters v characterizing the PSD. For a 
bimodal distribution, we also include the partial deriva- 
tive with respect to the fractional weight / in the second 
category. Thus, we distinguish two types of analytic 
derivatives: 

Type 1: with respect to single-particle characteristics: 


dij/ d\jj 
drrir ’ drrii 


(T-Matrix, Mie); 


di// 


(T-matrix only) 


Type 2: with respect to particle size distribution 
parameters and the fractional weight /: 

(T-matrix, Mie) 

dvk df 

Flere, is any one of the PSD parameters. For 

spheroids, shape factor £ is the ratio of the two semi-axes 
(oblate, £ > 1; prolate £ < 1; sphere £ = 1); for cylinders, £ 
is the diameter to height ratio; for Chebyshev particles, £ 
is the deformation parameter [29]. Some remarks are in 
order: 


(1) Mie scattering can be formulated as a special case of 
the T-matrix theory. It is possible with the NASA-GISS 
T-matrix code to obtain results for spherical particles 
to a high degree of accuracy by using a limiting 
case for spheroidal particles for which £ takes a value 
very close to 1.0 [9]. In practice, it is better to use a 
dedicated stand-alone Mie code for applications 
requiring spherical particle scattering, and there are 
a number of codes available in the literature. In this 
paper, we have created a stand-alone linearized Mie 
package in tandem with the linearized T-matrix 
model. 

(2) We do not consider derivatives with respect to the 
equivalent-sphere radius. For a single particle this 
radius is an input parameter; for polydisperse parti- 
cles, equivalent-sphere radii are specified through 
the PSD function. However, when the “equivalent- 
surface-area-sphere” representation is used in the 
T-matrix code, it is necessary to calculate the particle 
surface area and volume. Both these quantities are 
functions of the shape factor £, and their derivatives 
with respect to £ must be factored into the computa- 
tion of overall optical property derivatives d^/de. 
These additional derivatives are not required for the 
linearized “equivalent-volume-sphere” representation 
in the T-matrix code. 
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(3) For bimodal polydisperse applications, derivatives 
with respect to the number density fractional weight 
/are trivial; indeed from Eqs. (6a) and (6b) we find 

^Cext,sca _ ^(1) _f"(2) ^7;!^ 

ext, sea '^ext,sca> 


[cSy^>-cg>yf>]-yjC-cg>] 

df /cS+(i-/)C® 


(7b) 


2.2. The T-matrix ansatz and its linearization 

For electromagnetic scattering by an arbitrary fixed 
homogeneous object, expressions for the incident, inter- 
nal and scattered electric fields (Emc, ^par and Esca^ respec- 
tively) in terms of vector spherical wave functions 
and Nmn [8] are 

Umax ri 

Emc(J^) = [tlmnRS^mn{h^)^bmnRS^mn{h^y\f (8) 

n = 1 m = -n 
thnax ri 

Epar(R) = ^ ^ ^ ^ [Cmn^S^mnifOch^) dmn^S^mn03lc^^^y\f 

n = 1 m = -n 

( 9 ) 

Esca(R) = Wmn^rnn(l<^) + Qmn^rnn(hE()]- (10) 

n = 1 m = -n 

Here, R is the radius vector with origin inside the particle 
(which has circumscribed radius ro), k is the wave number 
and nic the complex refractive index of the particle 
(relative to the outside medium). Linearity of Maxwell’s 
theory and the boundary conditions dictates that there 
must be a linear relationship between the incident {amn, 
bmn] and scattered {pmn, Qmn] held coefficients; we express 
this in terms of the T-matrix T 


T-matrix formulation is really useful. The rotational 
transformation rule for the T-matrix is [30] 

"TLn..n.= E E 

mi = -nrri 2 = ~n' 

(14) 

Here, are the Wigner D functions, and (a, p, y) the 
Euler rotation angles. The pre-suffices on the T-matrix 
entries denote coordinate systems 1 and 2. This is an 
important result; once the T-matrix is known in coordi- 
nate system 1, then Eq. (14) allows us to calculate it in 
any other system. For rotationally symmetric particles, a 
convenient system takes the z-axis as that for rotation, 
and in this system the T-matrix has the symmetry relation 

'^mnrn'n' = ^mm'T^mnmn' 

We also note the relation T^Lm'n' 
which is a consequence of scattering matrix reciprocity 
[8]. For particles with spherical symmetry, the T-matrix 
ansatz reduces to 


T-l 1 ^ fi • t'22 

^ mnvn'n' ~nnn'i-fn, i mnm'n' 


—Snn' On, 


T-l 2 

^ mnm'n' 


t-21 

^ mnm'n' 


= 0 . 
(15) 


Here, a„ and are the usual Lorenz-Mie coefficients. 

Eq. (14) is the basis for averaging over particle orienta- 
tions. For the case of randomly oriented particles and the 
incident field in the form of a plane electromagnetic wave, 
the Wigner D-function orthogonality property allows us 
to derive the following well-known results for the extinc- 
tion and scattering cross-sections for randomly oriented 
particles [8]: 

rimax n 

Cext = --yReE E + (16) 

n = 1 m = -n 


Umax rimax ri n' 2 2 2 

E E EEh'-.-h 

n = \n' = im = -nm' = -n' i = \j =\ 


(17) 


p 

— T- 

a' 


'pi 1 

'pl2 


a 

q 


b 


'Y’21 

p-22 


b 


Similarly, one may write down linear systems relating 
the incident and internal fields, and the scattered and 
internal fields 


a 


Q" 

Q’2 


c 


P 


■RgQ" 

RgQ’2- 


c 

b 





d 


q 



Rg<i“ 


d 


(12) 


Combining (12) and (11), we find 


'pll 'pl2 


■RgQ" RgQ'2- 



1 

■t 

to 




q21 q22 


or (13) 

In this expression, matrices RgQ^ and Q are constructed 
from vector spherical wave functions that have been 
integrated over the particle’s surface. These spherical 
functions are products of well-known analytic functions, 
based on Bessel and Wigner d functions; detailed formu- 
lae are given below. 

Averaging over orientations is essential for non-sphe- 
rical particles, and it is here that the analytic nature of the 


For computing orientation averages of the scattering 
matrix expansion coefficients, the approach follows the 
use of the Clebsch-Gordan expansion for the Wigner d 
functions [8]. This has proved convenient for the compu- 
tation of nested sums of T-matrix coefficients; more 
details in Section 3.1. 

We now consider the T-matrix linearization. If x is 
any of the Type-1 variables (refractive index component, 
shape factor), then we may differentiate Eq. (13) directly 
to obtain the derivative T-matrix: 

ffT/dx = -5[RgQ,]/ax ■ Q,"’ -RgQ, ■ 5[Q."’ ]/5x (1 8) 

Now, since Q, Q,“’ =E (the identity matrix), we find 

a[Q-’]/5x = ■a[Q]/5x-(i-’ (19) 

Substituting (19) in (18), we find 

Sl/dX = -{d[Rg(iydX+l-d[(l]/dx) ■Q,”’. (20) 

The major computational task in determining the 
T matrix is evaluation of the inverse matrix in the 
NASA-GISS code, the LAPACK software is deployed for this 
task. We see in Eq. (20) that the only additional work 
required for computing the linearized T-matrix is the 
determination of derivatives of the matrices RgQ and Q, 
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since and T itself are already available to us. Compu- 
tation of these derivative matrices will also be dealt with 
in the next section. 


3. Type 1 derivatives for the T-matrix and Mie codes 

3.1. Vector spherical wave functions 


For the scattered field in Eq. (10), the vector spherical 
wave functions Mmn and are given by 

Mm„(kR) = (-l)'”d„h<„’>(x)C„.„(9)e™«’, x = kR; 

(21) 

N„„(kR) = hm(x)P^„(S) + 1 A [xf,m(x)] B„ 



(22) 

Bmn(^) ^ ^ ^ [^0m(’^)] + 9 ^Om(^)’ 

(23a) 

Cmn(^) = ^ dom(’9)-(p ^ [dom(^)]’ 

(23b) 


3.2. Linearization of Bessel functions 


Bessel functions in the Mie and T-matrix codes are 
determined by recursion. If x is the (real-valued) radial 
coordinate kR, then the downward recursion for spherical 
Bessel function j„(x) is given by 

F„(x) = Xj„(x); f n(x) = Gn(x)f „_1 (x), (29a) 


G„(x) = 


'2n + l 

X 


n -1 


~Gn + i(X) 


Gni (^) — 0- 


(29b) 


Here, Ni is the recursion starting value; there are a number 
of ways of setting this point. We use the specification in 
[4,8], namely, Ni(x) = x+4.05x^/^ + 60. For linearization, 
there is no dependence on refractive index variables, but 
X will depend on the particle shape parameter s if we are 
using the equivalent-surface-area-sphere (ESAS) represen- 
tation. Thus we must also consider the linearized recursion 

F'n(x) = xfnix) +xj;(x); f ^(x) = G„(x)F;_-, (X) + G;(x)Fn_i (X); 

(30a) 


G'„(x) = [G„(x)]2 


G'„+i(x)+ 


2n + l ; 


C'nS’O-O. 


(30b) 


P„„(a) - I dL(3); d„ - (23c) 

Hankel functions (of the first type) are h\l\x), and the 
Wigner d functions are given by 

(24) 


n (-1)"-"' / (n + m)! 

2" y (n-/)!(n + /)!(n-m)!' 


(25) 


Here, fi = cos9. Eq. (24) is valid for n>n* = max(|l|,|m|); 
otherwise djj^(9) = 0 for n < n*. The relations between the 
Wigner d and D functions and the generalized spherical 
functions are 

= i"-'PL(cos 9). 

(26) 


The orthogonality condition for the Wigner d functions is 
£ dl,^mt^(P)smpdp = ^ . (27) 


The orientation averaging proceeds through use of the 
Clebsch-Gordan expansion: 


d" TF)d"' .(B)= V 

^ ^nmn'xxix 


m + mi,m' + m'-^ 


iff 


(28) 


For details, see [8,27]. 

For the fields in Eqs. (8) and (9), the RgM and RgN 
functions are obtained by replacing the Hankel functions 
h^jl\x) by Bessel functions j„(x) and y„(x). For the interior 
field (Eq. (9)), we require Bessel functions of complex 
argument z = nicX = (irir + imOx. 


Prime indicates derivative d/ds. The determination of 
dx/ds is given in the next sub-section. 

Similarly the upward recursion for spherical Bessel 
function y„(x) is 

F„(x) = -xy„(x); F„+i(x)= +^F„(x)-F„_,(x); (31a) 

F_i(x) = sinx; Fq(x) = cosx. (31b) 

In this case, we use Ni as the recursion finishing point. 
The linearization with respect to the particle shape para- 
meter £ proceeds in a similar fashion 

F'„{x) = -x'y„(x)-xy'„(x); (32a) 

= ^^^^^lxF'„iX)-x!Fn(xy\-F'„_^(x); 

F'_^(x) = cosx; Fo(x) = -sinx. (32b) 

For complex-valued Bessel functions, we require a 
downward recursion similar to Eqs. (29a) and (29b), 
except in place of particle size parameter x, we have the 
complex argument z = (rur + m^x. For the refractive index 
linearizations, we then have 

F^(z) = z'Cniz) +zc;(z); F;(z) = Gn(z)F;_-i (z) + G;(z)Fn_i (z); 

(33a) 


G'^(Z) = [Gn(Z)f 


Gn+i(z) + 


2n + l 


G^2 = 0- 


(33b) 


Here, the prime symbol indicates derivatives d/drur or 
idldrup Similar considerations apply to d/d£. The recursion 
start is defined similarly through N 2 (z) = z+4.05z^/^ + 60. 

Mie formulae. In this special case, we can go directly 
to the Lorentz-Mie On and bn coefficients through the 
well-known results: 

_ [^ + ?]n(x)--F„-i(x) 

[^ + “]4>„(x)-<F„_i(x)’ 
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u [mcKn(Z)+g]n(X)-yn-l(X) , 

" [rricK„(z) + 2 ] (X) ’ 

= y„(x) = x/„(x); <f„(x) = x[;„(x)-iy„(x)]. 

(34b) 

Differentiations with respect to and rrii follow directly 
by chain-rule application of the above formula. 


3.3. Surface integral linearization (T-matrix only) 


Surface integrals are discussed in detail in [9]; here 
we summarize key formulas and focus on linearization 
aspects. For non-spherical particles the radius r = r(S,(p) is 
a function of angular coordinates S and (p, and we must 
therefore consider quantities such as 

jTn(r) ■ {RgMm'n’(kr,9,(p) x Mmn(kr.9,(p)]dS. (35) 


when evaluating the T-matrix. These surface integrals are 
calculated in spherical coordinates, with n(r) the outward 
normal vector. In general we may write for unit vectors 
{f.S.cp) [30]: 


n(r)dS = 


\ 1 dr - 1 dr / 

^ r dS rsinSdcp^ 


r^sin,9di9d(p. 


(36) 


In our case with rotationally symmetric particles, there 
is no azimuth dependence, so that r = r(d) only and the 
term for vanishes. Thus we need to evaluate the 
difference of two integrals: 


Spheroids. The explicit form used here is 


r(S) = 


V sin^i9-F£2 cos2,9 


1 dr _ (£2-i)sindcosd 
r dd sin^d-Fs^cos^d ’ 


( 39 ) 


Here, R(s) is the equivalent sphere radius, and shape factor 
£ is the ratio of the vertical and horizontal semi-major axes. 
For the equivalent volume sphere, R(s) is not dependent on 
£, but for the equivalent surface-area sphere, this depen- 
dency must be accounted for (Section 3.4 below). 

Linearization will require differentiation of Eqs. (38) 
and (39) with respect to the shape factor. Differentiating 
through the integrals in Eqs. (38a) and (38b) we find 



NG 


JJwin 

i = 1 




(40a) 


ds 



Wj 



Differentiation of Eq. (39) yields 

dr(S) ^ r 1 dR(s) ^ 1 £COS^d 1 ^ 
de R(s) d£ 3£ sin^S+s^cos^S 


(40b) 


(41a) 


d /I dr\ _ 2£sindcosd r41h^ 

d£ V ddy [sin^ cos^ 

As noted already, the derivative of ^(£) in Eq. (41a) will 
be zero for the equivalent volume sphere. 

Cylinders. There are two surfaces here, and the shape 
factor is now the diameter to height ratio. Then the 
quadrature is split according to [9]: 


p 71 p 71 

Jr= F(r,d)r^ sin ddd; Js= G(r,d) ^r sin ddd. 
Jo Jo 


(37) 


r(S) = (2/3)’/^6’/^4^; 

^ Sind 


1 dr(S) 
r dS 


-cotd (tand > £); 

(42a) 


Here f(r,d) and G(r,d) are the r and 6 components 
respectively of the kind of cross-product vector terms 
seen in equations of type (35); exact forms need not 
concern us here. As noted in [8,30], these integrals are 
restricted to ranges [0,ti/2] for particles with a plane of 
symmetry perpendicular to the rotation axis; for example 
r(7i-d) = r(d) for spheroids. The integrals in Eq. (37) are 
done using a double-range Gaussian quadrature scheme 
{-/rj,Wj} and {-F/ij,Wj} over the half-range intervals [-1, 0] 
and [0,1], respectively, where ft = cosd and i = 1,2. . .NG/2. 
Thus we may write 

NG 

Jr= (38a) 

7-1 i=l 




(38b) 


These integrals depend on the type of particle. The 
choice of NG is critical to the convergence of the T-matrix 
solution; an initial value is chosen such that NG = LN^ax, 
where L is an integer (dependent on particle choice) and 
Njnax is the size of the matrix Q, For more details, see [9]. 


r(9) = (2/3)G3a-2/3 = +tan9 (tanS < e). 

(42b) 

The linearization with respect to s is easy; the non- 
zero terms are 


8e ^ ’ R(e) de 3a 


(tanO > 6); 


(43a) 


dr(d)_ r 1 dm 2' 

ds ^ ^ R(s) ds 3S 


(tand < £). 


(43b) 


Chebyshev particles. These particles are generated 
through continuous deformation of a sphere of radius 
ro = R(s) using a Chebyshev polynomial of degree n; the 
deformation parameter s is always less than one: 


r(9) = R(a) [1+ 6 cos n3] ; 1 ^ 


snsinnS 


(44) 


1 -F£Cosnd' 

The linearization with respect to s is also straight 
forward 


-L£ cos nd] +R(s)cosnS; 
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d /I dr\ nsinnS 

ds\r dSj [1 + £ cos nSf 


(45) 


3.4. Equivalent surface area sphere (ESAS) linearization 
(T-matrix only) 


In this section, we look at the non-zero derivative 
dR(s)/ds which applies in the ESAS representation. In this 
case, the equivalent sphere radius Rq (a free parameter 
that does not depend on the nature of the particle under 
consideration) must be multiplied by a factor S(£) which is 
related to the particle surface area and volume. Specifi- 
cally 


R(s) = RqS(s) — Rq 


Ev{o) . 
Ea(s)’ 


Ev(s) = 


V(£) 


(4/3)7T_ 


1/3 


Ea(s) = 


A(sy 


4n 


1/2 


(46) 


The volume and area functions V(£) and A(s) are treated 
separately for the three particle types here. Note that 
S(£) = 1 for the sphere. 

Prolate spheroids (£ < 1 ). The volume is V(£) = 
(4/3)710^5; the surface area and function S(£) in (46) are 


A(£) = 2na^H(s) = 2na^ 


1 + 


sin ^ V 1-£2 
£\/1-£2 


S(£) = V2£-^3^(£)-^/2. 


(47) 


Here, we have polar and equatorial radii a and b, respec- 
tively, such that s = alb. The derivatives are 




= -S 


3e~^2H ’ 


dH _ -£\/ 1 -£2 -(1 -2£2)sin V1-£2 
ds £2(1— £2)3/2 


(48) 


Oblate spheroids (£> 1). The volume is again V(£)i 
(4/3)7ia^b. Functions A(£) and S(£) are given by 


A(£) = na^H(s) = na^ 


2 H- 


1 , £-F V£2-1 

= ln 


£V£2-1 £-V£2-1 


function S(£) and its derivative are given by 



as(£) 

ds 


= 


(£- 1 ) 

3s(23-s)’ 


(50) 


Chebyshev particles. Radius r(S) is given by Eq. (44). To 
find the surface area and volume, we use a quadrature: 

NS I 

A(£) = ^ Wj[l -h£COsn,9i]Y [1 s cos nSif 4-^ s^n^ sin^ nSii 

i = 1 

(51a) 


NS 

V(£) = ^ WjSini9i[l 4- s cos nSif (sin di[\ 4- s cos nSi] 4- snxi). 
i= 1 

(51b) 

Here, the quadrature is {x,, W/} over the interval [-1, 1], 
and Xi = cosdi. The number of quadrature points NS =60 
in the original F77 NASA-GISS code [9]; we have retained 
this number. Differentiation of Eqs. (51a) and (51b) with 
respect to the deformation parameter is a lengthy but 
straightforward exercise. 


4. Type 2 derivatives for the T-matrix and Mie codes 


For polydisperse applications, we use a range of 
(equivalent-sphere) particle size distributions present in 
the Mie code of [4]; details, see Table 5. Suppose i/^(r) is 
any monodisperse optical property to be integrated over 
size radius. Then the PSD integrations are done using a 
series of Gauss-Legendre quadratures {rkj,Wkj}J< = 
1, . . .NQO), one for each block j, where there are NB blocks 
covering the full range [ri,r 2 ] 


1 EffiEKn(rfej.v)^(rkj)Wkj 


(52) 


If Vq is one of the set v of (up to 3) parameters character- 
izing the PSD, then the linearization of Eq. (52) with 
respect to Vq is 


dVq 


^NB 
2^j = l 


dn(ruj,v) 
2^k = 1 dVq 
si^NB s^NQij) 
= 1 2^k = 1 


n(rkj,v)wi,j 


(53) 


S(e) = (49) 

Linearization with respect to s proceeds by analytic 
differentiation of Eq. (49); this is a straightforward alge- 
braic exercise. 

Cylinders. Here, factor s is the diameter to height ratio. 
The formulas for this case are particularly simple; the 


It is well known that orientation averaging in the 
T-matrix solutions tends to reduce high-frequency varia- 
tions, so that it is not necessary to use more than one 
quadrature block in the PSD integration. Thus for T-matrix 
polydispersion, NB=1. 

One example will suffice to illustrate the PSD linear- 
ization process. For the lognormal distribution with 


Table 1 

Sample finite difference validation. 



Analytic T-matrix Jacobians 


Einite-difference T-matrix Jacobians 



•2<^(Cext) 

2{(Csc) 


5^(Cext) 

g{(C5ca) 



1.0228E + 01 
1.0343E-02 
-4.8533E-01 

8.9886E + 00 
-6.4536E-01 
-2.8879E-01 

-1.6865E-01 

3.7953E-02 

-7.7285E-02 

1.0226E+01 

1.0343E-02 

-4.8522E-01 

8.9872E + 00 
-6.4536E-01 
-2.8869E-01 

-1.6857E-01 

3.7953E-02 

-7.7280E-02 
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parameters Tg (mode radius) and Sg (standard deviation) 
we have 


n(r) = 


V^rsc 


exp 


(Inr-lnrg)^ 


2s| 


(54a) 


an(r) 

drg 


m 

h 


(Inr-lnrg); 


5n(r) 

dSg 


n(r) 


1 - 


(\nr-lnrgf 


(54b) 


Similar results can be established for the gamma, mod- 
ified-gamma, and power-law distributions commonly 
found in the literature; in all cases, analytic differentiation 
of these well known functions is straightforward. 


5. Some results 

5.1. Finite difference testing 

All linearized outputs from the monomodal T-matrix 
code are normalized, that is, if we are seeking a Jacobian 
with respect to quantity c^, then the actual output is (say 
for the extinction coefficient) 

£^(Cext)^^^. (55) 

Given this definition, one way to obtain a finite difference 
estimate of the derivative is 

rr \ ^ ^ext(i )—Cext(i) 

0(^((-ext) = g * (5b) 

where the perturbed quantity is ^ = ^{\+5) for some 
small number 6. Comparing (55) and (56) allows us to 
make a finite-difference validation of the analytic weight- 
ing functions in a convenient manner. In practice the 
optimum value of 6 will depend on the parameter under 
consideration; experience with this testing indicates that 
d = \0~^ is best for the shape factor and refractive index 
real part parameters, whereas (5 = 10“^ is good enough for 
the refractive index imaginary part and the PSD para- 
meter derivatives. 

Table 1 has examples of finite difference Jacobian 
validations for the extinction and scattering cross-sec- 
tions and the asymmetry parameter, with calculations for 
monodisperse oblate spheroids with mr+im/= 1.42 + 
0.005Z, shape factor s = \.l, and particle size parameter 
12.56; a perturbation 6 = \0~^ gives results accurate to 
the 4th significant figure. 

For a single-mode call to the linearized Mie code, all 
weighting function outputs are unnormalized (absolute 
derivatives). For the bimodal applications, T-matrix deri- 
vatives are normalized and Mie derivatives again unnor- 
malized. The only exception is with fractional number 
density derivatives d<i/^>/d/, which are unnormalized for 
both codes. 

Finally we note that the computer code package has a 
facility for carrying out this finite difference validation for 
any type of particle, as well as some coding for testing the 
new T-matrix Fortran 90 against the old NASA-GISS 
Fortran 77 package. 


5.2. Examples of output 

Here we present some sample results, focusing on the 
linearized optical properties. This section is intended to 
give a flavor of the kind of output generated by the 
linearized model; specific retrieval applications are beyond 
the scope of the present work. In order to give an overview, 
we have employed color contour plots similar to those for 
example in [5] (Plates 2.1 -2.4). We look at the following 
situations for randomly oriented rotationally symmetric 
particles. 

In Fig. 1, we look at the extinction cross-section Qxt 
and its two normalized derivatives mrdCextIdmr and sdCextl 
ds for oblate and prolate spheroid particles, with incident 
light at wavelength 0.55 pm and fixed imaginary refrac- 
tive index component mi =0.005. Results are plotted for 
a range of values [0.4, 2.0] for the shape factor s, and 
a range [1.1, 1.6] for the real refractive index component 
m^. Calculations were done using the equivalent surface 
area sphere (ESAS) representation. Left panels show results 
for a monodisperse situation with particle size 1 pm (size 
parameter ~ 12.56), with the right panels containing 
results for a polydisperse aggregate characterized by a 
lognormal PSD with mode radius 0.5 pm and standard 
deviation 2 pm. 

Focusing next on Chebyshev particles in Fig. 2, we look 
at the extinction cross-section Cext and the single scatter- 
ing albedo co (top left and top right, respectively), and 
their normalized derivatives midCextIdmi and midcoldmi 
(middle row) with respect to the imaginary component 
mi of the refractive index, and derivatives sdCextIds 
and sdcolds (bottom row) with respect to the Chebyshev 
deformation parameter. Incident light has wavelength 
0.95 pm and the real part of the refractive index compo- 
nent is fixed at mr= 1.33. Results are plotted for a range of 
values [0.01, 0.3] for the deformation parameter £, and a 
range [0.002, 0.22] for the imaginary refractive index 
component m/. Calculations were done using the equiva- 
lent surface-area-sphere (ESAS) representation. 

In Fig. 3, we return to oblate spheroids, looking this 
time at angular distributions. Results are shown for 
monodisperse spheroids with shape factor £ = 1.7 and 
refractive index 1.42 + 0.008/, at wavelength 0.443 pm. 
We look at the normalized scattering matrix element 
fii(0) and the corresponding degree of linear polariza- 
tion (in %) -F 2 i( 0 )/Fii( 0 ) (top left and top right, respec- 
tively), along with their three derivatives sd/ds, m^jdmr 
and midjdmi (rows 2 to 4, respectively). Results are plotted 
against scattering angle 0 from 0 to 180°, and for a range 
[0, 20] for the particle size parameter. Calculations are 
again done using the ESAS representation. The plot for 
-F 2 i( 0 )/Fii( 0 ) (top right) is closely similar to one of the 
graphs in Plate 2.1 of [5]. It can be seen that polarization 
at backscattering angles has large sensitivity to changes of 
shape factor; this offers the promise for retrieval of 
particle shape from multi-angle polarization measure- 
ments. The sensitivity drops signiflcantly at scattering 
angle close to 180° for particles with size parameters less 
than 5; this is also the condition for low polarization. 

In Fig. 4, we look at a bimodal aggregate, comprising 
a dust mode with lognormal polydisperse spheroidal 
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Spheroids C_ext, Size parameter 12.56, nj 0.005, Monodisperse(left), Lognormal PSD (right) 
Normalized Jacobians: W.R.T. n_r (middle row), W.R.T. EPS (Bottom row) 



Real (Refractive Index), n_r 


Real (Refractive Index), n_r 


Fig. 1. Spheroid extinction cross-sections Cext (top panels) plus derivatives nirdCextlOnir (middle panels) and sdCextiOs (lower panels), contour-plotted 
against shape factor e and real refractive index component nir. (Left column) monodisperse particles with size parameter 12.56; (right column) 
polydisperse lognormal aggregates with mode radius 0.5 pm and standard deviation 2 pm. Fixed imaginary refractive index component m, = 0.005, 
equivalent surface-area-sphere (ESAS) representation. 


particles with mc= 1.53 -fO.OOSi, and a sulfate mode with 
lognormal polydisperse spheres of refractive index rric= 
1.43 -F 2 X 10“^i. We look at the scattering matrix element 
Fii (plotted on a natural logarithmic scale for conveni- 
ence) and the degree of linear polarization Pun (in %) 
-F2i(0)/fii(0) (first and second rows), along with their 
sensitivities with respect to the fractional number density 
weight /of the dust mode (rows three and four). Results 
are plotted against the fraction of dust for a fixed dust- 
particle shape factor 1.7 (left column), and against the 
dust particle shape factor (varying from 0.7 to 2) with a 
fixed fractional weight of f=0.5 (right column). The 
sensitivities are here defined as the (normalized) deriva- 
tives of the Ln(Fii) or Pun with respect to/. Variations of 
both Fii and Pun with respect to/(i.e. vertical changes of 
color in the upper two panels on the left) are much less 
than their counterparts with respect to shape factor 
(i.e. color changes in the two upper panels, right column). 
An overview of the sensitivity of such variations with / 
and shape factor is seen in the lower two rows of Fig. 4; 
moreover. Pun has relatively larger sensitivity to / than 
Fii. Note that, since Pun is less than 40% (second row in 
Fig. 4), the normalized relative sensitivity /d[Ln(PnNj]/d/ is 
larger than fdPuNidf as shown here. [fd[Ln{PuN)]/^f is not 
shown here, as Pun is zero at scattering angles 0° and 180]. 
Overall, Fig. 4 suggests that, at least for the aerosol para- 
meters specified here, angular polarization is useful for 
retrieving the fraction of non-spherical large particles [31]. 


6. Computer codes 

The initial-release package of linearized FORTRAN 90 
T-matrix and Mie codes may be obtained upon inquiry 
from the corresponding author; the codes are in the 
public domain, and when the codes become optimized 
and better established, it is intended that they will again 
be available from the GISS website (http://www.giss.nasa. 
gov/~crmim). The F90 T-matrix package is based on the 
existing NASA-GISS F77 code [9], while the Mie package is 
based on the Meerhoff code [4]. The new F90 codes are 
accompanied by a User Guide. 

For the T-matrix part of the package, the following 
remarks apply to the F90 upgrade: 

• Most of the original naming conventions in the F77 
code have been preserved; 

• all subroutines have declared input and output explicitly 
— no common block storage; 

• all code is “implicit none” with explicit declaration of 
all variables; 

• all subroutine I/O has explicit intent (In/Out/lnOut) 
signifiers; 

• equivalence statements have been removed; 

• code has an explicit exception handling procedure for 
dealing with program failure; 

• dimensioning is symbolic throughout, no allocatable 
arrays (at least in this version). 
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Monodisperse Chebyshevs (T2), size par = 6.61, n_r 1.33, C_ext (left), SS_alb (right) 
Normalized Jacobians: W.R.T. nj (middle row), W.R.T. EPS (Bottom row) 
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Fig. 2. Extinction cross-sections Cext (top left) and single scattering albedos co (top right) for monodisperse Chebyshev particles with fixed real refractive 
index component mr= 1.33 and size parameter 6.61. ESAS calculations were done for a range of deformation parameters [0.01, 0.3] and a range of values 
of the imaginary component of the refractive index from 0.002 to 0.22. Normalized derivatives niidCextlOmi and m,5co/5mi are shown in the middle panels, 
with associated (normalized) deformation parameter derivatives sdCextiOs and sdcolds in the lower panels. 


The choice of PSD for the F90 T-matrix code has been 
extended to include all the options present in the Meerhoff 
[4] code (see Table 5). However, all original PSD specifica- 
tions from the F77 NASA-GISS code have been preserved, 
and the user can still choose one of these PSDs by turning 
on the Boolean flag “Do_PSD_oidstyle” (see below). This 
option is useful for validation against the F77 code. Note 
that the PSD-linearization is only possible with the new- 
style PSD choices, for which the PSDs were explicitly 
differentiated as part of the linearized Mie package. 

The code has been made more flexible, with a greater 
range of input choices, now specified by reading from 
configuration files — an example of this is described in detail 
below for the T-matrix case. As far as output is concerned, 
some users may require just the “bulk” optical properties 
(extinction and scattering cross-sections, single scattering 
albedo), and in this case, the code calculating the expansion 
coefficients is turned off. Note however that if the asym- 
metry parameter is desired, then the expansion coefficient 
code must be activated. Similarly, in many applications (e.g. 
when providing optical property inputs for radiative transfer 
modeling) it is only necessary to compute the bulk quan- 
tities and expansion coefficients — thus the F-matrix output 
for a regular grid of scattering angles is optional. 

6A. Code descriptions 

There are two main directories in the package (Fig. 5). 
The main Tmatrix_environment directory contains a 


number of “makefiles”, which will generate executables 
named according to any test programs present. Modules 
and object files are stored in separate subdirectories to 
avoid clutter. This directory also contains the configura- 
tion files. Results files may also be stored separately. The 
Mie part of the package is structured similarly. 

The complete set of modules to be used in any call 
to the linearized T-matrix model is found in the other 
directory Tmatrix_sourcecode, which contains the 14 files 
outlined in Table 2. The Mie code is simpler; all functions 
(with the exception of the parameters, the bimodal 
masters and the I/O read and write routine) are contained 
in the two master modules (Table 3). 

6.2. Configuration file example 

Table 4 contains an example of a configuration file. 
Each file is divided into 4 groups. In this particular case, 
the program will perform a full calculation of bulk 
properties and expansion coefficients and the F-matrix, 
using the equivalent surface area sphere (ESAS) represen- 
tation with a lognormal PSD (new-style) of mode radius 
0.5 pm and standard deviation 2 pm, with limiting radii 
0.1 and 5 pm, for oblate spheroids of shape factor 1.7 and 
refractive index (1.53, 0.005) at wavelength 0.55 pm. This 
input will also generate 5 linearizations: weighting func- 
tions with respect to the two refractive index compo- 
nents, the shape factor and the two PSD parameters. 
Similar configuration files have been designed for the 
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Left Column : Ln(F_11) ; n_r.d(Ln(F11 ))/dn_r; nJ.d(Ln(F11))/dn_i; eps.d(Ln(F11))/deps 
Right Column: PJin = -F_21/F_11 (in %); n_r.d(P_lin)/dn_r; n_i.d(P_lin)/dnJ: eps.d(P_lin)/deps 



scattering angle (degs) scattering angle (degs) 

Fig. 3. Logarithm of the (1, 1) scattering matrix element Ln(Fn(0)) (top left panel); degree of linear polarization (in %) -F 2 i( 0 )/F„( 0 ) (top right panel). 
The three derivatives nirdldnir, midldnii and sdids of these quantities are shown in rows 2, 3 and 4, respectively. Results are plotted for monodisperse 
oblate spheroids (e = 1.7, mr= 1.42, m, = 0.008) against scattering angle 0 for a range of particle size parameters as indicated. Calculations were again 
done in the ESAS representation. 


Mie code. In addition, there are configuration files for 
bimodal applications (requiring two sets of microphysical 
values and PSD inputs). 

6.3. Exception handling 

Both the T-matrix and Mie codes have consistent 
exception handling procedures for dealing with input 
checking and execution failures. An overall Boolean flag 
is output for “success/failure”, and there is also an integer 
status variable plus three character strings for output 
messages. Inputs are checked for consistency, and in the 
case of input error, a message will be generated describ- 
ing the error, plus a second message outlining the action 
required to correct the error along with 2 or 3 traces to 
establish the location of the error. Dimensioning and 
convergence issues are the main causes of T-matrix 
execution failure, and the appropriate messages from 
the F77 code have been retained in the F90 package. In 
both codes, dimensioning checks will suggest new para- 
meters to use. 

6.4. Particle-size distributions 

Table 5 summarizes PSD options in the T-matrix and 
Mie packages. For PSD_Index=3 (Old-style power law) 
and FixRlR2=T, then PSD_Parl is an effective radius, 
PSD_Par2 an effective variance. 


7. Concluding remarlcs 

In this paper we have described a complete lineariza- 
tion of the T-matrix model as it applies to randomly 
oriented axially symmetric particles. The linearity of 
Maxwell’s equations and the intrinsic analytical nature 
of the T-matrix code allow us to carry out analytic 
differentiation of the entire T-matrix solution with 
respect to any variables characterizing the particles in 
question. We distinguish two types of linearization: (1) 
with respect to single particle characteristics (real and 
imaginary components of the refractive index, particle 
shape or deformation factor), and (2) with respect to 
particle size distribution parameters characterizing poly- 
disperse aggregations. 

The NASA-GISS T-matrix code package has been trans- 
lated to Fortran 90, and additional code written to gen- 
erate the Type 1 and Type 2 optical property derivatives 
as noted in Section 2.1. The new code has been validated 
against the old NASA-GISS FORTRAN 77 package, and all 
optical property derivatives have been checked against 
flnite-difference estimations. We have developed a sepa- 
rate linearization package for the Mie code, even though 
Mie theory is a special case of the T-matrix formulation. 
The Type-1 Mie linearization applies only to derivatives 
with respect to the refractive index components. Type-2 
linearizations apply equally to the Mie and T-matrix 
formulations. 
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Fig. 4. (First row) Logarithm of the (1, 1) element Ln(fn(0). (Second row) Degree of linear polarization (in %) Pun^ -f2i(0)/fii(0)- (Third and fourth 
rows) Derivatives /5[Ln(Fii)]/5/ and /5 [Pujv]/ 5/, where / denotes the dust aerosol volume fraction. Results are for bimodal aerosol particles with a 
lognormal spherical sulfate mode (mr=1.43, m, = 2 x 10“^, rg=0.04, cTg=1.8) and a lognormal spheroid dust mode (mr=1.53, m, = 0.005, rg=0.4, cTg=1.8). 
Results are plotted against scattering angle 0 for a range of dust-mode aerosol number density fraction (from 0.01 to 0.96) with fixed shape factor 
(e = 1.7) in the left column, and for a range of shape factors (from 0.7 to 2) with fixed dust-mode volume fraction (/=0.5) in the right column. Calculations 
were again done in the ESAS representation. 



Fig. 5. Directory structure of the linearized T-matrix/Mie package. 
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Table 2 

Subdirectory “Tmatrix_sourcecode”, F90 modules. 


Module 


Purpose 


tmat_parameters . f 90 
tmat_distributions . f 90 
tmat_master . f 90 
tmat_master_PLUS . f90 
tmat_f unctions . f90 
tmat_functions_PLUS . f90 
tmat_makers . f90 
tmat_makers_PLUS . f90 
tmat_scattering . f 90 
tmat_scattering_PLUS . f 90 
Utilities_LAPACK. f90 
Tmat_IO_readwrite . f 90 
Tmat_master_bimodal . f90 
Tmat_master_bimodal_PLUS . f90 

Dimensioning and type-kind parameters 

Computation of PSDs (old style and new style) and linearizations (new style only) 
Top-level module for computing optical property stuff, standard output only 
Top-level module for computing optical property standard and linearized output 
Work-horse routines for Bessel and other functions 

Work-horse routines for Bessel and other functions, and any linearizations thereof 
Routines for creating Q, RgQ, and T-matrix 

Routines for creating Q, RgQ. and T-matrix, and all linearizations thereof 
Routines for calculating expansion coefficients and F-matrices 
Routines for expansion coefficients and F-matrices, and all linearizations thereof 
LAPACK routines (F90 syntactical translation of original F77 code) 

Routines for reading configuration files and writing standard and extended outputs to files 

Bimodal wrapper, standard output 

Bimodal wrapper, standard and linearized output 

Table 3 

Subdirectory “Mie_sourcecode”, F90 modules. 

Module 

Purpose 

Mie_parameters . f90 
Mie_distribution . f90 
Mie_main . f 9 0 
Mie_main_PLUS . f 90 
Mie_IO_readwrite . f 90 
Mie_master_bimodal . f90 
Mie_master_bimodal_PLUS . f 90 

Mie dimensioning and type-kind parameters 
PSD distribution functions 

Mie module for computing optical properties, standard output only 

Mie module for computing optical properties, standard linearized output 

Routines for reading configuration files and writing standard and extended outputs to files 

Bimodal wrapper, standard output 

Bimodal wrapper, standard and linearized output 


Table 4 

Configuration file example for the T-matrix model. 


Value 

Name 

Description 

Remarks 

*** First group (Boolean flags) 


T 

Do_Expcoef f s 

Flag for expansion coefficient output 

New feature 

T 

Do_Fmatrix 

Flag for optional F-matrix output 

New: Do_Expcoeffs must be set 

F 

Do Monodisperse 

Flag for a monodisperse calculation 

New feature 

T 

Do_EqSaSphere 

Flag for using equivalent surface area sphere (ESAS) 
representation 

Formerly a non-Boolean input 

T 

Do_LinearRef 

Flag for linearizing w.r.t. real and imaginary parts of 
refractive index 


T 

Do_LinearEps 

Flag for linearizing w.r.t. shape parameter 


T 

Do_LinearPSD 

Flag for linearizing w.r.t. PSD parameters 

Only works for the “New-style” PSD choices 

F 

Do_psd_OldStyle 

Flag for using original PSD choices 

If set, use the NASA-GISS F77 original PSD choices 

*** Second group (PSD control) 


4 

psd_ index 

Particle size distribution (PSD) index 

See Table 5 for choices 

0.5 

psd parsl 

First PSD parameter 

See Table 5 

2.0 

psd pars2 

Second PSD parameter 


0.0 

psd_pars3 

Third PSD parameter 


1.0 

Monoradius 

Size (pm) of equivalent-sphere particle 

Monodisperse only 

F 

FixRlR2 

Flag for fixing R1 & R2 internally 

Only if Do_psd_01dStyle not set 

0.1 

R1 

Minimum radius (microns) 

Not needed if FixRlR2 is set 

1.0 

R2 

Maximum radius (microns) 

Not needed if FixRlR2 is set 

*** Third group (General control) 


-1 

np 

-1 (spheroids), -2(cylinder), > O(Chebyshev) 

Same as G1SS-F77 options 

20 

nkmax 

Number of PSD quadrature points 

Same name as in G1SS-F77 

91 

npna 

Number of F-matrix outputs 

Same name as in G1SS-F77 

2 

ndgs 

Number of ESAS division points 

Same name as in G1SS-F77 

2.0 

eps 

Aspect ratio, deformation parameter , etc. 

Shape parameter 

0.001 

accuracy 

Accuracy for convergence 

As in G1SS-F77, formerly DELT 

*** Fourth group (optical inputs) 


0.5 

lambda 

Wavelength 

Always micrometers 

1.53 

n_real 

Real part of refractive index 


0.008 

n_imag 

Imaginary part of refractive index 
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Table 5 

Summary of PSD options. 


Old Tmatrix New Tmatrix/Mie Description of PSD_lndex 


4 


2 


1 Two parameter GAMMA with ALPHA and B given 

2 Two parameter GAMMA with REFF and VEFF given 

3 3-parameter bimodal equal-weight GAMMAS, 2 REFF, 1 VEFF 

4 Two parameter LogN with RG and SIGMA given 

5 Two parameter LogN with REFF and VEFF given 


3/5 6 

7 


Power-law with Rl, R2 and ALPHA (=3, old style) 
3-parameter modified-gamma with ALPHA,GC,GAMMA given 


1 


8 


3-parameter modified-gamma with ALPHA,B,GAMMA given 


At present, the code is restricted to double precision 
floating point arithmetic; in the next version we plan to 
allow for a user-specifled level of numerical precision 
(this is a nice F90 feature), so that the code can be run in 
“extended precision” mode without the need for a sepa- 
rate package as is currently the case with the GISS F77 
code. This will extend the usage to particle size para- 
meters in excess of 100. Performance and allocatable- 
memory optimizations are also planned. 

All software in the package is in the public domain; 
the codes may be downloaded on a trial basis from RT 
Solutions by contacting the corresponding author. This is 
the first “beta” version; user feedback will help to con- 
solidate this code and improve portability and robustness, 
and it is intended that the second release will be made 
from the NASA GISS website. 

The original motivation for the development of this 
package has come from Earth-atmosphere remote- 
sensing inverse problems for aerosol retrieval, with the 
emphasis on retrieving microphysical aerosol character- 
istics rather than macrophysical optical properties. In this 
context, the package is best used in conjunction with a 
linearized radiative transfer model such as VLIDORT; such 
a combination is then able to deliver forward-model 
analytic Jacobians necessary for aerosol retrieval pro- 
blems using least-squares fitting (with or without reg- 
ularization). T-matrix codes for non-spherical scattering 
are found in many other atmospheric physics applications 
as well as diverse fields such as hydrometeor scattering 
and biomedical applications, and it is hoped that the 
present code will contain something for everyone. 

We note that there is no reason why the linearization 
process described here cannot be extended to other non- 
spherical scattering situations using the T-matrix approach 
(for instance, coated or chiral particles, situations with 
non-random orientations, etc.). Future work will focus on 
the linearization of these variants of T-matrix theory, and 
we will also focus on the generation of Jacobians from 
distributions of homogeneous spheres. 
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